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Abstract 

Effective modeling and numerical spectral-based propagation schemes are proposed for ad- 
dressing the challenges in time-dependent quantum simulations of systems ranging from atoms, 
molecules, and nanostructures to emerging nanoelectronic devices. While time-dependent Hamil- 
tonian problems can be formally solved by propagating the solutions along tiny simulation time 
steps, a direct numerical treatment is often considered too computationally demanding. In this 
paper, however, we propose to go beyond these limitations by introducing high-performance nu- 
merical propagation schemes to compute the solution of the time-ordered evolution operator. In 
addition to the direct Hamiltonian diagonalizations that can be efficiently performed using the new 
eigenvalue solver FEAST, we have designed a Gaussian propagation scheme and a basis transformed 
propagation scheme (BTPS) which allow to reduce considerably the simulation times needed by 
time intervals. It is outlined that BTPS offers the best computational efficiency allowing new 
perspectives in time-dependent simulations. Finally, these numerical schemes are applied to study 
the AC response of a (5,5) carbon nanotube within a 3D real-space mesh framework. 
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I. INTRODUCTION 



Nowadays, the numerical solution of the time-dependent Schrodinger-type equation is still 
considered as one of the most challenging problems in quantum simulations of molecules, 
nanostructures and devices. In nanoelectronic applications, in particular, efficient time de- 
pendent simulations have become increasingly important for characterizing the electron dy- 
namics under time dependent external perturbations such as electromagnetic fields, pulsed 
lasers, AC signals, particle scattering, etc. For example, THz frequency responses for carbon 
nanotube (CNT) have recently been observed in [lj pointing out the need of time-dependent 
simulations capable of going beyond the linear response regime. Reliable modeling ap- 
proaches in time domain, however, are often limited in term of trade-off between robustness 
and performances 13] . 

One approach for solving the time-dependent Schrodinger-type equation consists of using 
a partial differential equation representation where one can generally discretize the time 
domain using finite difference method. The specific techniques include both explicit and 
implicit schemes, with the commonly used Crank-Nicolson scheme [4j. These numerical 
techniques can be cast as direct approaches, however, they can end up being numerically 
expensive in the case of long time quantum simulations where it is essential to preserve 
accuracy and robustness. 

Another commonly used approach consists of solving the integral form of the problem 
via the numerical treatment of the time-ordered evolution operator. Two cases can then 
be generally considered: (i) The Hamiltonian is time-independent; (ii) The Hamiltonian is 
time-dependent. In order to solve this latter, one traditional solution consists of dividing 
the simulation time domain into tiny time steps A t: such that the Hamiltonian can be 
considered as time-independent within A t and the system can be solved with techniques 
similar to case (i). Case (i) is indeed formally straightforward, since the problem is then 
equivalent to solving the exponential of a Hamiltonian. The most obvious way to address 
this numerical problem would be to directly diagonalize the Hamiltonian while selecting 
the relevant number of modes (i.e. eigenpairs) needed to accurately expand the solutions. 
Spectral decomposition, however, are known to be computationally demanding especially for 
large systems. The mainstream in time dependent simulations has then been to avoid solving 
directly the eigenvalue problem and use approximations most often based on split operator 
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techniques [5HZJ - In this paper, however, it is pointed out that direct diagonalizations can 
now be efficiently performed by taking advantage of the recent high-performance eigenvalue 
solver FEAST [SI Eg. 

However, when the Hamiltonian is time-dependent, a direct propagation technique still 
involves solving a very large number of eigenvalue problems along the whole simulation times. 
Although, the FEAST solver can be used at each time step to speed up the numerical process, 
we propose here to introduce novel spectral-based propagation schemes that can significantly 
reduce the total number of eigenvalue problems while preserving high numerical accuracy 
and robustness. The high-efficiency of these techniques can be achieved, in particular, by 
considering larger A t and focusing only on obtaining the final states at the end of each time 
interval. 

This paper is organized as follows: Section [TT] presents the basics of time dependent 
Schrodinger equation within the TDDFT and Kohn-Sham formalism; Section [TTT] describes 
the various propagation approach used in this work including a direct scheme, a novel Gaus- 
sian quadrature scheme, and an optimized basis transformed propagation scheme (BTPS); 



The numerical efficiency of these techniques is finally compared in Section IV using simula- 
tion results on 3D CNT. 



II. TIME-DEPENDENCE IN QUANTUM SYSTEMS 

In quantum systems, the electrons obey the time-dependent Schrodinger equation: 

ih^m = mt). (i) 

Using a single electron picture, the Hamiltonian is composed of two terms, one of kinetic 
origin and another describing the interaction of the particle with a local potential which 
may be time dependent: 

H = f + V(t) = -^-V 2 + v(r,t). (2) 

Besides appropriate boundary conditions, the time dependent Schrodinger equation requires 
an initial value condition ty(t = 0) = that completely determines the dynamics of 
the system. In particular, within the time dependent density functional theory (TDDFT) 
framework [10J, the solutions of the stationary Kohn-Sham Schrodinger-type equation are 
taken as initial wave functions and will be propagated over time. 
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ih—ipj(r, t) 



h 2 

-—V 2 + v KS [n](r,t) 
2m 



V>;M). (3) 



TDDFT can be viewed as a reformulation of time dependent quantum mechanics where 
the basic variable is no longer the many body wave function, but the time dependent electron 
density n(r, t). For a system composed of N e electrons, the electron density can be obtained 
from the solution of a set of one body equations, the so called Kohn Sham equations, that 
have the same form as Equation ([I]) where ^ = {^i, . . . , ^N e } and with ipj solution of: 

d_ 

The density of the interacting system can be obtained from the time-dependent Kohn-Sham 
wave functions 

N e 

n(r,t) = ^|^-(r,t)| 2 . (4) 

The Kohn-Sham potential Vks is a functional of the time-dependent density and it is con- 
ventionally separated in the following way: 

v KS [n](r,t) = v ext (r,t) +v H [n](r,t) + v xc [n](r,i), (5) 

where the first term represents the external potential, the second term is the Hartree po- 
tential which accounts for the electrostatic interaction between the electrons, and the last 
term is defined as the exchange-correlation potential which accounts for all the non-trivial 
many-body effects. Finally, it should be noted that at t = 0, the initial wave functions 
= ^2°\ • • • > ^iv e } are solutions of the ground state DFT Kohn-Sham stationary 

equations [IT] : 

~ + VKs[n}(v)] < } (r) = Ef^°\r). (6) 

In a confined system, the Kohn-Sham wave function ip^ (i.e. ^-(r, £ = 0)) is associated 
with the eigenvalue Ej°\ and the ground state many-body density is given by n(r, t = 0). 
Formally, the solution of Equation ([!]) can be written as 



= U(t, 0)^o = T exp |~ J\tH{t)^ * , 



(7) 



where the evolution operator U is unitary and can be represented using a time ordered 
exponential Texp, which is a non-trivial mathematical object defined in noncommutative 
algebras. It is important to note that if the Hamiltonian is time-independent, the solution 
takes a simplified form: 

^(t) = expj-^#W (8) 



Unfortunately, this is not the case in most relevant applications which require the description 
of the electron dynamics under time dependent external perturbations such as electromag- 
netic fields, pulsed lasers, AC signals, particle scattering, etc. In this case, approximations 
such as perturbation theory or linear response are commonly used to simplify the compu- 
tational cost of the time-dependent solutions. In the following Section |III[ however, novel 
spectral-based propagation schemes are proposed to address the numerical challenges of 
solving ([7]) for the case of large scale problems and long simulation time domains. In order 
to ease the description of the numerical techniques, one will consider only non-interacting 
systems (i.e. the potential v is time dependent but it is not a functional of n). The natural 
extension of these numerical schemes to solving non-linear problems will be discussed in 
Section [V] 



III. SPECTRAL-BASED PROPAGATION SCHEMES 
A. Direct propagation approach 

In practice, intermediate physical solutions are computed in addition to the final solution 
in order to describe the evolution of the system over [0, £]. This can be accomplished by 
dividing [0, t] into smaller time intervals since using the intrinsic properties of the evolution 
operator, one can apply the following decomposition: 

U(t, 0) = U (t ni * n -i)l7(*n-l, *n- 2 ) ...U (f 2 , *l) , * ) , (9) 

where we consider n—1 intermediate interval time steps with £ = 0, and t n = t. Most often 
a constant time step A t is used and one has to deal now with the problem of performing 
many shorter time propagation of the solutions along the whole interval [0,£]: 

^(t + A t ) = Texp | — ^ jf + * drH(T)}v(t). (10) 

Additionally, if A t is chosen very small, it is reasonable to consider H{r) constant within 
the time interval [t,t + A t ] leading to 

(f + A,) = U(t + A t , t)*(t) = exp \- l -A t H(t)\ (11) 



h 

which is also equivalent to solving the time independent problem ([8]). 



Denoting H the N x N Hamiltonian matrix obtained after the discretization of H at a 
given time t and where TV could represent the number of basis functions (or number of nodes 
using real-space mesh techniques), H can then be diagonalized as follows: 

D = P T HP, (12) 

where the columns of the matrix P = {pi,P2, • • • ,Pm} represent the eigenvectors of H 
associated with the M lowest eigenvalues regrouped within the diagonal matrix D = 
{di, d 2 , . . . , dyi}- If the discretization is performed using non-orthogonal basis functions 
(e.g. finite element basis functions in real-space), the eigenvalue problem that needs to be 
solved takes the generalized form: 

H Pi = diSpi, (13) 

where S is a symmetric positive-definite matrix. Here, we consider that the computed 
eigenvectors P are S-orthonormal i.e. P T SP = I. In order to evaluate the exponential of 



the Hamiltonian in (11), it is now possible to perform a spectral decomposition where the 
exponential acts only on the eigenvalue matrix D. One can show that the resulting matrix 
form of time propagation equation is given by: 



*(t + A t ) = P exp |-^A t D j P T S*(t), (14) 

which is exact if M = N. We also note that * T S* = I with * = {^i, V>2j • • • , ^ Ne }- I n 
practice, it is reasonable to obtain very accurate spectral approximations even if one selects 
the number of lowest eigenvalues M much smaller than the size of the system TV but greater 
than N e (i.e. N e < M « N). The N e propagating states * are indeed low energy states. 
In this case, P, D, S are respectively matrices of size TV x N e) N x M, M x M, and 
N x N . Using the property ([9]), the solution \l/(£) can finally be obtained in function of 



Pi exp ( -^A t Di)PfS 



*o, (15) 



*(*)=r|n 

where = Pf H(tj)Pj, and the time ordering is defined by 

T jn A{U) | = A(t N ) . . . A{t 2 )A(t x ). (16) 

Therefore, the direct approach may require solving hundred to thousand of eigenvalue 
problems all along the time domain (one eigenvalue problem for each time step). For large- 
scale systems (e.g. where hundreds of atoms are taken into consideration), this approach 
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is often considered not applicable as accurate eigensolutions for thousands of eigenpairs are 
computationally demanding and algorithmically challenging. Recently, however, the new 
density-matrix based algorithm FEAST [8] has been proposed to overcome these difficulties, 
and provide both accuracy, robustness and performance scalability for solving eigenvalue 
problems. FEAST is a general purpose algorithm which takes its inspiration from the 
density-matrix representation and contour integration technique in quantum mechanics. Its 
main computational tasks consist of solving very few independent linear systems (typically 
eight linear systems that can be solved simultaneously in parallel) and one reduced dense 
eigenvalue problem orders of magnitude smaller than the original one (of size MxM, in the 
present case). FEAST is also ideally suited for the direct propagation scheme presented here, 
since it can take advantage of the subspace computed at a given time step i as initial guess 
for the next time step i + 1 in order to speed-up the numerical convergence. Ultimately, one 
can show that if enough parallel computing power is available at hand, the main computa- 
tional cost of FEAST for solving the eigenvalue problem can be reduced to solving only one 
linear system. In contrast to the proposed FEAST approach for the direct spectral-based 
propagation scheme, PDE-based numerical schemes for solving ([7]) (such as Crank-Nicolson) 
would require solving successively a large number of linear systems by tiny time intervals 
A t . A high-performance implementation of the FEAST algorithm can be found in [9]. 



B. Gaussian propagation scheme 

The direct approach has a sounded physical interpretation as it corresponds to a step by 
step propagation of the solution over time. Its major drawback, however, is that it involves 
a very large number of time steps from initial to final simulation times. Let us now consider 
the case of a much greater time interval A t which may correspond, for instance, to a given 
period of an external time-dependent perturbation. The solution can then be obtained by 



performing first a numerical integration for the integral on the Hamiltonian in (10): 



+ A t ) = Texp j-l^^tf^) j (17) 

where uji and £ are integration constant, and p is the number of quadrature points. If the 
quadrature points are close enough, one can additionally assume that the anti-commutation 
error between Hamiltonians can be ignored, and the exponential can then be decomposed 



into a product of exponentials for each time step 



. 2=1 



*(*)• 



(18) 



The number of exponentials to evaluate along [0, t] is then equal to the number of time 
interval n multiply by the number of quadrature points p by intervals. After discretization 



and spectral decomposition of each Hamiltonian as presented in Section III A, it comes: 



*(i + A t ) = r n 



. i=l 



P i exp(--^Di )PfS 



*(*)■ 



(19) 



We note from the derivation of equations (17) and (18) that two numerical errors are 
respectively involved: (i) an integration error and (ii) an error on the anti-commutation 
resulting from the decomposition of the exponential. The direct propagation scheme pre- 
sented in Section |III A| can be derived by assuming a rectangle quadrature rule with uJi — 1 
and using £ = 8 t , U = t + (i — 1) * St where 5 t = A t /(p — 1) corresponds to a very small 
time step within intervals. Using this less conventional derivation, it is possible to point 



out that a drastic reduction of the integration error (17) could be obtained by considering 
higher-order quadrature scheme such as Gaussian quadrature. A p-point Gaussian quadra- 
ture rule is a numerical integration constructed to yield an exact result for polynomials of 
degree 2p — 1 by a suitable choice of the points U and (Gauss-Legendre) weights u^. We 



associate the quadrature points ti at the Gauss node X{ using ti 



+ ^±At, also we 



note £ = At/2. Gaussian quadrature can use relatively much fewer points than low-order 
quadrature rule such as rectangle, to yield a high order approximation of the integral of a 
function. Interestingly, a large number of numerical experiments have shown that accurate 



exponential decomposition (18) can still be obtained while using much fewer points p or 
larger spacing within intervals (even in the presence of strong perturbations). The same 
experiments have demonstrated, however, that increasing the accuracy on the numerical 
integration plays an important role for obtaining the correct solutions. We will show in 
Section [IV| that the Gaussian quadrature scheme provides a good combination of reduced 
computational consumption and high numerical accuracy. The number of exponentials that 
needs to be evaluated by time intervals can indeed be reduced by a factor 3 to 4 as compared 
to the direct approach. Using the Gaussian quadrature, however, the intermediate solutions 
obtained by propagating the solutions within the intervals have no physical meaning. Since 
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the quadrature weights uji are different for each Gauss node, the solutions can indeed only 
be known at the beginning and at the end of each interval. 



C. Basis transform propagation scheme (BTPS) 

By focusing on obtaining only the final states at the end of each time interval A t , the 
Gaussian propagation schemes can considerably reduce the number of eigenvalue problems 
that needs to be solved. Yet, solving large scale eigenvalue problems is still the most compu- 
tational consuming part in our simulation. Here, we introduce a basis transform propagation 
scheme (BTPS) which can help in reducing not only the number of eigenvalue problems to 
solve by intervals, but also the size of each eigenvalue problem. 

Before performing the decomposition of the exponential, let us first consider the dis- 



cretization of equation (17): 

*(t + A t ) = Texp j-l^^H; J (20) 

where = H(^). We propose then to project the "pseudo-Hamiltonians" cc^H^ onto a 
common eigen-subspace P constructed from the result of the quadrature sum ^f =1 ^H^. 
Denoting H this new "global Hamiltonian" which takes the overall contribution of the time- 
dependent perturbation over a given time interval A t (however H has no physical meaning) , 
one can perform the diagonalization defined in (fl2j) i.e. 



D = P T HP = P T UiH-i \ P, (21) 



i=i 



where the subspace P and eigenvalues D are obtained solving the eigenvalue problem in 



(13). Denoting hi the projection of the pseudo-Hamiltonians at a given quadrature point 
such that 

h, = P T (^H,)P, (22) 



it comes with (21) 



D = J>,, (23) 



1=1 



where are M x M dense matrices while their summation D is diagonal. Using a spectral 
decomposition on the time-ordered exponential for the global Hamiltonian, the expression 
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(20) can now take the following form: 

*(* + A t ) = PT exp \ -k h i \ P T S*(t), 



(24) 



2=1 



where P has been moved outside the expression of the time-ordered exponential since it is 
a common subspace for all different quadrature time points within [t,t + A t ]. It is then 
important to note that one cannot replace directly the sum of the projected Hamiltonian 



in the exponential by D (23), as the time-ordered operator can only be resolved using a 



product of functions taken at different times (16). Similarly to the derivation of (18), we 



assume that the quadrature points are close enough that one can decompose the exponential 
into a product of exponentials for each time step: 

p 



¥(f + At) = PT \ JJexp \ --£h, |> \ P T S*(t), 



.2=1 



(25) 



Thereafter, one can perform the diagonalization of each hj matrices 



(26) 



where the M x M matrices represent the eigenvectors of h^, and the diagonal matrices 
ei regroup the associated M eigenvalues. In contrast to the generalized eigenvalue problem 
on the global Hamiltonian above, these eigenvalue problems are standard i.e.: hq = qe. 



Finally using a spectral decomposition on the exponentials in (|25|), we obtain: 

p 



*(t+A ( )=pr m 



.i=i 



qiexp 



h 



&i qi 



P T S*(t). 



(27) 



In contrast to the propagation schemes presented in Sections HI A| and |III B which required 
solving n *p large scale eigenvalue problems of size N x N along [0, t] (n number of intervals 
At, and p number of quadrature points), the BTPS approach consists of solving only one 



NxN eigenvalue problem by time intervals (21 ), and n*p reduced dense eigenvalue problems 



of size M x M (26). 



The values of £, U and uji have been specified in Section |IIIB| in the case of a rectangle 
and Gaussian quadrature rules. We recall that the Gaussian scheme has been proposed to 
reduce the number of quadrature points p (i.e. number of exponential to evaluate) within 



interval. From equation (24), however, one can see that the numerical integration is here 



actually performed on the reduced projected Hamiltonians, and this basis transformation 
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is expected to significantly decrease the integration error. As a result, we have found that 
the rectangle quadrature rule is likely to provide the same accuracy than the Gaussian 
scheme using the same number of quadrature points p. Since this latter cannot be reduced 
below a certain threshold in order to limit the error arising from the decomposition on the 



exponentials (25), a high-order integration scheme such as Gauss quadrature becomes then 
obsolete using BTPS. Similarly the Gaussian scheme, however, intermediate solutions that 
can be obtained by introducing I = P T SP within the product expression in (27), would 
end up having no physical meaning. Indeed, one can demonstrate that using BTPS the 
(S)-orthonormalization of the wave functions is only satisfied between the solutions taken 
at the beginning and at the end of each interval. 



IV. SIMULATION RESULTS 



We propose to perform time-dependent 3D simulations of an isolated carbon nanotube 
(CNT) which is sandwiched between two electrodes producing a AC voltage at the THz 
frequency. In this model, the charge transfer at the contacts is not considered (i.e. open 
systems and transport problems are not considered). Here, our model uses a local empiri- 
cal pseudopotential approach, real-space mesh techniques for discretization (finite element 
method), and a time-dependent version of the atomistic mode approach described in | [T2l [T3] . 
Moreover, the empirical pseudopotential U eps is supposed time-independent the total atom- 
istic potential can then be decomposed as follows: 

U(x,y,z,t) = U eps (x,y,z) + U ext (x) sin(cjt), (28a) 

2r — L 

U ext (x) = ~^Uo, (28b) 

where x is the longitudinal direction of the tube, L represents the distance between contacts 
(x e [0,1/]), oj = 27r/, with / the corresponding frequency of the AC signal and Uo its 
amplitude. The time-dependent external potential applied to the CNT, maintains then zero 
in the middle of the tube but oscillates at both ends alternatively to the ±Uq values. 

Here, the system under consideration is a 6-unit cell of a (5,5) CNT where the contact- 
contact distance is set at 1.98nra. This CNT is composed by 120 atoms, and using the 
empirical pseudopotential for Carbon atom proposed in [14], one can expect to capture all 
the 60 electrons (without spin-dependence) contained in the tt z orbitals. Using the notations 
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defined in Section III, it comes N e = 60, and M = 120 (M being the number of eigenpairs 
used for the spectral decomposition and chosen here as twice the number of electrons N e ). 
The finite element discretization of the full 3D system gives rise to sparse matrices of size 
N ~ 500* 10 3 , while this size can be drastically reduced in our case using a mode approach in 
real-space to N ~ 15* 10 3 . For the external perturbation (AC signal), we consider Uq = 5eV 
and / = 200THz. In our simulation, the solution wave functions \l/ = . . . , i/) Ne } will 

be propagated from t = to t = 8T where T — 1/ f — 5 * 10 -15 s denotes the period of the 
AC signal. The time interval used is A t = T where we have found that p = 120 integration 
points by intervals are sufficient for the direct propagation scheme to provide an accurate 
reference solution. Let us recall that direct scheme allows to capture the evolution of all the 
wave functions at each tiny time step A t /p, while the Gaussian and BTPS approaches are 
expected to produce only accurate results at each time interval A t which length has been 
arbitrarily chosen very large here to point out the robustness of the approaches. 

In order to examine the relative accuracy and robustness of the numerical techniques 
presented in this article, we propose to calculate the energy evolution of the wave functions 
i.e. 

Ej(t) = <^,(t)|H(t)|^,(t)>, j = l,...,N e . (29) 



Let us first discuss the integration error which is introduced in our model in (17). As an 
example, Figure [T] compares the direct and Gaussian propagation schemes by representing 
the energy evolution of both the first energy level E\ and the HOMO level of the system Eqq. 
Using a direct scheme with p = 40 integration points for performing the rectangle quadrature 
between time intervals and as compared to the reference curves obtained for p = 120, one 
can see that the results start diverging after few time steps. In contrast to the direct 
scheme, accurate results can still be obtained at the end of each time period using p = 40 
and higher order integration technique such as our proposed Gaussian scheme. In practice, 



Gauss quadrature could accurately perform the numerical integration (17) using much fewer 



integration points (p < 40), however, this would also increase the anti-commutation error 



arising from the decomposition of the exponential in (18) as discussed in Section III B and 
shown in Figure [2} 

Finally, similarly to the Gaussian propagation scheme, p = 40 integration points suffice 
for the BTPS approach to accurately obtain all the solutions at the end of each time period. 
Using BTPS the numerical integration is performed on a projected space, and this operation 
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FIG. 1: Evolution of the energy expectation of the first level (on the left) and the HOMO level (on 
the right) along 8 time periods of the AC signal. The straight lines represent the reference solutions 
calculated using the direct scheme with p = 120 rectangle quadrature points within intervals. The 
results for the energy evolution using a direct scheme p = 40, represented using dashed lines, 
diverges after few time steps, while the same number of points did suffice for the Gaussian scheme 
to accurately capture the solutions at the end of each time period (represented by filled triangle 
symbols). 



does not require high-order integration schemes as discussed in Section [Til C[ As compared to 



direct or Gaussian propagation schemes, BTPS reduces drastically the computational costs 
meanwhile stable and accurate. However and as shown in Figure |3j intermediate solutions 
that can be computed within the intervals have no physical meaning. It should be noted 
that the same remark applied to the Gaussian scheme, although the (unphysical) solutions 
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Time (5*10 s) 



FIG. 2: Evolution of the energy expectation the levels 57 to 60. The filled symbols are obtained 
using p = 40 Gaussian propagation scheme, and match exactly the reference solutions (dashed 
lines) at the end of each time period. In contrast, the solutions obtained using the p = 16 Gaussian 
scheme, represented by the unfilled symbols, suffer now from an unsuitable approximation on the 



decomposition of the exponential (18) due to an increase in distance between integration points. 



within intervals were not represented in Figures [T] and [2] for clarity. 

From the solutions on the wave functions, one can now investigate many properties of 
the CNT. In particular, within the real-space mesh framework, the electron density is given 
by: 

N e 

12 (30) 



As an example, the results of the 3D simulations on the electron density are represented in 
Figures [5] and |5] respectively at t = and t = T/A using a contour plot. 
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FIG. 3: Evolution of the energy expectation for the HOMO level. BTPS can capture efficiently and 
exactly the reference solution at the end of each time interval using p = 40 points for the rectangle 
quadrature rule (the BTPS solutions are pointed out in the plot using the circle symbols). We 
note that the intermediate solutions obtained with BTPS have no physical meaning. 

In order to better illustrate the time evolution of the density, it is possible to calculate the 
variation of the ID projection of the electron density on the longitudinal axis i.e. rimfe) = 
J y z n(x, y, z, t)dydz. Figure |6] shows how this ID electron density calculated using the direct 
propagation scheme and at some particular positions along x, evolves over time. It should 
be then noted that much fewer points that the p = 120 required by A t using the direct 
approach, could obviously be sufficient to accurately capture the variation on the density. 
Indeed, the electron density as well as other integrated quantities, are likely to exhibit much 
weaker variations as compared to the variations of the individual wave functions. Let us then 
suppose that the variation of the electron density can be captured using (at most) 30 points 
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FIG. 4: Contour plot showing the 3D electron density at t = 0. 

by A t = T, BTPS could then be efficiently used to compute these intermediate solutions 
using a new time interval of A t = T/30. For the same degree of accuracy, indeed, the direct 
propagation scheme would now require solving 4 (120/30) large eigenvalue problems within 
the new A t , but BTPS would still require solving only one (and 40/30 c± 1 small reduced 
eigenvalue problem). 

V. CONCLUSION 

In this work, we have investigated three different spectral-based propagation techniques 
for time dependent quantum system: Direct, Gaussian and BTPS approaches. These nu- 
merical schemes have been applied to study the AC response of an isolated single wall 
carbon nanotube using a real-space mesh techniques framework, empirical pseudopotential 
and non-interacting TDDFT calculations. Using the direct approach, the time-ordered evo- 
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FIG. 5: Contour plot showing the 3D electron density at t = T/4, electrons have moved to the left 
edge. 

lution operator is solved via a step-by-step diagonalization procedure of the time-dependent 
Hamiltonian. Spectral-based schemes are known as robust and accurate but are traditionally 
considered too computationally expensive for addressing large systems. The new eigenvalue 
solver FEAST, however, can efficiently address these problems and provide performances 
and scalability. We have also pointed out that two numerical errors do appear in time prop- 
agation problems: a quadrature error resulting from the integral on the Hamiltonian, and 
an error resulting from the decomposition of the exponential. In contrast to the more con- 
ventional direct approach, our proposed Gaussian scheme demonstrates that it is possible to 
obtain the solutions at the end of each time interval using a reduced number of quadrature 
points (i.e. reduction of the quadrature error using a higher-order integration scheme), while 
preserving accurate exponential decompositions. Finally the BTPS scheme reduces not only 
the number of eigenvalue problems to solve by intervals, but also the size of each eigenvalue 
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3.5 




FIG. 6: Evolution of the ID projection of the electron density at two given positions. In this 
example, the positions xiand X2 are selected between contacts and edges of the CNT, and the 
results are then representative of the CNT field emission. 

problem. Since only one large scale eigenvalue problem and a small number of reduced 
eigenvalue problems need to be solved by time intervals, the computational efficiency and 
time savings offer by BTPS are significant. 

It is straightforward to note that the direct scheme can also be used to address self- 
consistent TDDFT calculations using the adiabatic local density approximation (ALDA) 
where the potentials vh and v xc ^ need to be calculated at each time step in function 
of the local density. In order to take advantage of the BTPS scheme for interacting sys- 
tems, however, one would need to adequately keep track of the variation of the density 
self-consistently with the local potential by considering a small enough time interval A t . 
Since density and potential are likely to exhibit much weaker variations with time as com- 
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pared to individual wave functions, BTPS should still provide significant computational time 
savings as compared to the direct approach. 

To summarize, the spectral-based propagation schemes proposed here with in partic- 
ular the optimized BTPS approach, are potentially capable to open new perspectives in 
time dependent simulations of large-scale quantum systems. Possible applications of these 
techniques range from obtaining accurately the excited states of arbitrary molecules and 
nanostructures, to efficient characterization of high frequency responses of emerging nano- 
electronic materials and devices. 
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